Spatial vector solitons in nonlinear photonic crystal fibers 
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• We study spatial vector solitons in a photonic crystal fiber (PCF) made of a material with the 

' focusing Kerr nonlinearity. We show that such two-component localized nonlinear waves consist 

, of two mutually trapped components confined by the PCF linear and the self-induced nonlinear 

refractive indices, and they bifurcate from the corresponding scalar solitons. We demonstrate that, 
' in a sharp contrast with an entirely homogeneous nonlinear Kerr medium where both scalar and 

vector spatial solitons are unstable and may collapse, the periodic structure of PCF can stabilize 
the otherwise unstable two-dimensional spatial optical solitons. We apply the matrix criterion for 
' stability of these two-parameter solitons, and verify it by direct numerical simulations. 
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PACS numbers: 42.65.Tg, 42.65.Jx, 42.70.Qs 



I. INTRODUCTION 



Photonic crystal fibers (PCF) have attracted much interest due to their intriguing properties, many potential 
applications, as well as the recent development of successful technologies for their fabrication with engineered linear 
and nonlinear properties P, Q . Photonic crystal fibers are characterized by a conventional cylindric geometry with a 
J> ' two-dimensional lattice of air-holes running parallel the fiber optical axis. Such PCF structures share the propagation 
properties of photonic crystals, based on the existence of the frequency gap with the transmission suppressed due to 
CO ' the Bragg scattering, as well as the properties of conventional optical fibers, due to the presence of a defect in the 
\ structure acting as a PCF core. Some of the PCF intriguing characteristics include the possibility to design single- 
moded PCFs independently on the light frequency even for a large core, allowing the guidance of high powers what 
makes PCFs very suitable for amplifiers or laser cavity applications. On the other hand, there exists an upper cutoff 
frequency by means of a reduction of the core index, and this also allows a very flexible control on the dispersion 
properties, supporting large shifts of the zero-dispersion point, and birefringence, which can be made much higher 
than in conventional fibers by a proper design. 

In PCFs, light confinement is restricted to the core of the fiber and therefore nonlinear effects, such as light self- 
trapping and localization in the form of spatial optical solitons 0, may become important. The stabilizing effect 
• of periodic media for optical solitons has been observed in a number of cases. In particular, one-dimensional vector 
solitons that are unstable in uniform media were made stable in a medium with a periodic modulation of the refractive 
index (jl. Also, discrete vector solitons where experimentally observed in two-dimensional optically-induced photonic 
d ' lattices Similar to the case of two-dimensional nonlinear photonic crystals j^, it has been recently demonstrated 
numerically that a PCF can support and stabilize both fundamental and vortex spatial optical solitons 0,0. In a 
sharp contrast with an entirely homogeneous nonlinear Kerr medium where spatial solitons are unstable and may 
collapse, it was shown that the periodic structure of PCF can stabilize the otherwise unstable two-dimensional spatial 
optical solitons. 

In this paper,^e make a further step forward in the study of nonlinear effects in PCFs, in comparison with the 
recent analysis and analyze the existence and stability of spatial vector solitons in PCFs. In general, vector 

solitons are defined as two-component mutually trapped localized beams whose properties may differ substantially 



from the properties of one-component scalar solitons [3|- In addition, two-dimensional vector solitons are known to 
be unstable in the nonlinear Kerr medium (see, e.g., Ref. In contrast, as we show in this paper, the periodic 
modulation of the refractive index in PCF provides an effective physical mechanism to stabilize the otherwise unstable 
two-dimensional spatial optical solitons. We study the stability of these two-parameter solitons and apply the matrix 
stability criterion that is then verified by direct numerical simulations. 

The structure of this paper is the following. First, In Sec. II we introduce our physical model that is characterized 
by an effective potential created by the PCF environment and also describe the nonlinear interaction between the 
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beam components. Then, in Sec. Ill we introduce our numerical method to find the classes of spatially localized modes 
existing in the nonlinear core of the PCF. In Section IV we describe the family of two-component spatial solitons. 
Finally, in Sec. V the stability of both one- and two-component solitons is analyzed. 



II. MODEL 



We consider a simple model of PCF that describes, at a given frequency, the spatial distribution of light in a 
nonlinear dielectric material with a triangular lattice of air holes in a circular geometry. We assume that the PCF 
material possesses a nonlinear Kerr response, and the hole at the center is filled by the same material creating a 
nonlinear defect, as shown in Fig. ^a,b). In the substrate material of the fiber, the linear refractive index is n^, 
whereas inside the holes it is n^. Air holes have radius r. We consider the case when the PCF core guides two modes 
or two orthogonal polarizations. In the nonlinear regime, the mutual interaction between these two modes is described 
by the system of coupled equations, 

+ AiPi + uaipi + V{x, y){5 + + Af|V'2p)V'i = 0, 
+ Ail^2 + na^2 + V{x, v){5 + \ij2? + Ai|^lP)V'2 = 0, 

oz 

where "01 and "02 are two components (or two polarizations) of the electric field, A is a transversal Laplace operator in 
{x, y), S ^ Ug — ria, and V{x, y) is an effective potential describing the defect and the lattice of holes in the transverse 
plane (x, y). We normalize V = 1 m the material outside the holes, and F = in the holes. The nonlinear incoherent 
interaction between the components is described by the parameter ^. 

To find stationary two-dimensional nonlinear modes of PCF, we look for the solutions in the form 

■01 (a;, y, z) = u{x, y)e''^\ ip2{x, y, z) = v{x, y)e'^'' , 

and obtain the following coupled system of z-independent differential equations: 

(iu — Am -I- UaU + V{x, y) {5 + + fiv^) u, 

(2) 

— Av + UaV + V{x, y) {5 + + M^^) v. 

The model ^ describes the stationary distribution of a two-component field in an inhomogeneous nonlinear medium, 
in a planar geometry. Without the external potential, the vector solitons in both one- and two-dimensional cases 
have been studied earlier However, the lattice of air-holes and the central defect break the radial symmetry of the 
problem, and the corresponding vector solitons are not radially symmetric. 
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III. NUMERICAL METHOD 



In order to find the solutions for nonlinear localized modes, we consider a rectangular domain of the {x,y) and 
apply a finite-difference scheme, taking n and m uniformly distributed samples of the variables x and j/, respectively, 
in order to cover all the domain. Denoting those samples as < i < n, and yj,l < j < m, at each mapped 
point {xi,yj) of the domain we consider the corresponding samples for all the functions defined in the equations: 
Uij = u{xi, yj) and, similarly, the second component Vij, and the potential Vij. Substituting these re-defined variables 
into the model and imposing homogeneous boundary conditions in all four edges of the domain, we obtain an 
algebraic nonlinear problem of 2 x n x m equations with the same number of unknowns and Vij. 

In order to make the notation more compact, the samples corresponding to different functions, which constitute 
n X m matrices, are rearranged concatenating the columns of the matrices to produce big column vectors of N rows 
{N — n X m), u, v, and V. Besides, we compact the vectors corresponding to both field components in a unique 
field vector, by concatenating one after another: q = (u-^|v-^)^, with 2N components qk- In that way, the algebraic 
nonlinear system can be written as Aq = 0, where A is the 2N x 2N matrix which depends on the unknown vector 
through the nonlinear terms, and we denote as A[q], so that the system of equations takes the form, 

A[q]q = 0, (3) 

being the rows of the matrix product 

= E,(A[q])fc,q„ (4) 
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so that the system is written as = 0, k = 1,2, . . . , 2N. The matrix A, even being huge in size, is in practice 
very sparse, and it differs from zero at the main diagonal, two diagonals next to the main one, and two more at the 
distance n from the main one (this four diagonals appear due to the coupling terms in the derivatives of the Laplace 
operator), and also two more diagonals at a distance N from the main one, due to the coupling between both field 
components. 

The nonlinear system of equations JSJ can be solved using the standard globally convergent Newton method 
which builds the solution iteratively from an initial guess q'^ in the form q' — q'~^ + Sq, I = 1,2,..., where the 
calculation of the so-called Newton step 6q at each iteration involves the solution of the linear system: 

J{6q) = -E, (5) 

where E is the vector obtained by substituting the last iterate into Eq. Q), and J is the Jacobian matrix defined 
as Jij — dEi/dqj, and also evaluated substituting the last known iterate. The Jacobian matrix presents a similar 
sparse structure as the matrix A, and it can be calculated analytically. Obviously, due to a huge size of the matrix 
J, the system can only be solved iteratively. Taking into account that for our particular problem the matrix J is 
symmetric, though in general indefinite, the SYMMLQ method ^2 proved to be successful. 

Some improvements would be possible in the method, taking the advantage of the system symmetries. In fact, 
due to the hexagonal lattice of holes, the field should be invariant under the rotation by the angles ^(7r/3), where I 
is an integer. It would make possible to solve the problem only in a circular sector of the amplitude 7r/3, imposing 
periodic boundary conditions at the borders and homogeneous in the radial direction. The number of points could be 
reduced in that case by the factor of six. Another approach, that also takes an advantage of the lattice periodicity, 
was developed by Ferrando et al. . 

IV. STATIONARY SOLUTIONS 

The presence of the external linear potential given by the central defect and the lattice of air-holes makes the system 
non-scalable and its radial symmetry broken. Therefore, the study has to be carried out by numerical methods. We 
solve Eq. ||2Jl numerically for both scalar (when one of the components vanishes, i.e. v — Q) and vector (or two- 
component) spatial solitons and obtain the stationary states of the nonlinear system. 

A. Scalar solitons 

For the scalar case, we assume that one of the components is absent (e.g., u = 0) and we study a single nonlinear 
equation of the nonlinear eigenvalue problem We find a family of the spatially localized modes-the so-called PCF 
spatial solitons-as a function of the mode propagation number (3. These results are similar to those earlier reported 
by Ferrando et al. , and the solution can be envisaged as the fundamental mode of the effective fiber generated by 
the combined effect of the PCF refractive index and the nonlinear index induced by the solution amplitude itself. 

Figures nja,b) show two examples of stationary, spatially localized solutions of the nonlinear model ||2Jl at w = 0, 
which describe scalar spatial optical solitons as nonlinear modes of PCF. The whole family of such one-parameter 
solutions can be characterized by the power P = J v?dxdy, that is plotted in Fig. ^c), where the points A and B 
correspond to the examples (a,b), respectively. The material parameters for the PCF are taken as = 1, r^s = 4 
and r = 0.75. 

First, we notice that these stationary solutions for scalar spatial solitons in PCF have been found earlier by Ferrando 
et al. Q, who also mentioned, without a proof, that such nonlinear modes are stabilized by the lattices of PCF holes. 
Indeed, it is well known that in the nonlinear focusing Kerr media without a nonlinearity saturation, the self-trapped 
optical beams are always unstable Q. This instability can manifest itself as the beam spreading, when the input 
power is lower than that of the soliton, or the beam collapse, when the power is larger than the soliton power. As 
has been mentioned earlier by Ferrando et al. such a soliton instability can be suppressed by the presence of the 
lattice of holes, because the external potential stops the beam spreading, as it happens in a conventional optical fiber, 
leading to the existence of a family of stable stationary beams. 

In order to demonstrate this feature, we follow the standard analysis of the soliton stability I3| and plot in Fig.^c) 
the soliton power as a function of the soliton propagation constant. A positive slope of this dependence indicates 
the soliton stability, as will be demonstrated below. In Fig. |21we present some related numerical simulations of the 
dynamics of a perturbed scalar soliton. Some of the stationary states are scaled by factors slightly higher and lower 
than the unity respectively, so as to induce an initial perturbation, and then propagated using a standard beam- 
propagation-method algorithm. The result is that the soliton behaves stably if its power remains below the maximal 
limiting power on Fig.^c). 
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FIG. 1: (a,b) Examples of stationary solutions for scalar PCF spatial solitons for /3 = 3.2 and /3 = 6, and (c) power dependence 
of the soliton family. Points A and B in the power diagram (c) correspond to the examples (a,b), respectively. Inset: the 
corresponding one-dimensional profiles at y = for both examples. 
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FIG. 2: Results of numerical simulations of the soliton dynamics in the scalar case. The initial stationary solution is perturbed 
by two amplitude scalings (indicated in the graphs): one is higher and the other one is lower than the unity. Shown is the 
maximum soliton intensity vs. the propagation distance. 



When the scaling factor is taken higher than unity, a stable propagation is observed for the solitons of low enough 
power, as seen in Fig.l^Ja). Nevertheless, for higher values of the power the soliton may collapse if the scaling factor 
is too large [Fig.[2Ib)], but it remains stable for a smaller scaling, see Fig.|2Ic). Further increase in the initial power 
results in collapse of the beam for any scaling factor (Fig. I^Jd)). 

When the scaling factor is taken smaller than unity, the lattice of holes stops the soliton spreading in all cases, so 
that the soliton propagates stable, as is illustrated in all cases presented in Fig. |21 
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FIG. 3: Examples of the two-component stationary solutions of the model Q for vector solitons in PCFs. Two cases correspond 
to the points A (/3 = 5, 7 = 4) and B (/3 = 5, 7 = 8) in the existence domain shown in Fig. 2] for fi = 2. In each column both 
components are shown in the plane {x,y), together with a one- dimensional s-cutoff profile to compare the field amplitudes. 




FIG. 4: Existence domain for the vector solitons in PCF in the plane (/?, 7), shown at two values of the coupling parameter 
fi. Labeled points correspond to two particular examples shown in Fig.|3 



B. Vector solitons 



Vector solitons in the coupled problem Q depend on both propagation constants (/?, 7) as well as the material 
parameters {6,fi). Some examples of the vector solitons in PCF are presented in Fig. |31 corresponding to the points 
A and B marked on the existence domain shown in Fig. 0] This domain, whose symmetry respect both parameters 
P and 7 is evident from the symmetry of both the equations of the model lO, is plotted in the plane (/3,7). The 
existence domain is limited by two lines at whose points (bifurcation points) the vector solitons originate from the 
scalar solitons; such curves can be regarded as bifurcation curves. When /i < 1, close to the lower bifurcation curve, the 
second component decreases becoming a linear guided mode of the soliton mode in the first (self-guided) component; 
the opposite case occurs close to the upper bifurcation curve where the role of the components is reversed. When 
/i > 1, we have the opposed situation with respect to the lower and upper bifurcation curves. The presence of an 
effective waveguide associated with one missing hole in the lattice is the reason that the propagation constants take 
a value different from zero, when the power vanishes; this threshold value corresponds to the eigenvalue of the linear 
mode guided by this effective waveguide in the lattice of air-holes. 

Similar to the scalar case, the presence of a periodic lattice of holes suggests that the vector solitons may become 
stable in this system. In this case, the vectorial nature of the system plays an important role to determine the portion 
of the domain where the solutions are stable. As follows from the next section, by applying the generalized matrix 
stability criterion, it is possible to determine the boundary between the stable and unstable regions. According to 
that, this boundary is the set of points that fulfill the marginal stability condition det(£') — 0, where Dij = dPi/d(3j 
while (Pi, P2) and (/3i, /32) are respectively powers and propagation constants for the components (m, v). In Fig.0|this 
boundary, as well as both regions of stability and instability are represented. A number of numerical simulations were 
carried out to test the stability of the solutions in each region. A standard beam propagation algorithm was used and 
the stationary solutions of the system were rescaled by a constant slightly higher that one, so that the peak amplitude 
of the fields initially raises over the peak amplitude of the exact stationary solution. For fields in the stability region, 
in spite of the initial growth of amplitude, it becomes stable after certain propagation distance. 
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V. SOLITON STABILITY 



Stability of scalar and vector solitons in the coupled NLS equation can be studied with the matrix stability criterion 
|l3l Il4j . Applications of the matrix criterion depend on the exact count on the number of eigenvalues of the matrix 
Schrodinger operators and require careful numerical computations of a spectral (linearization) problem. Alternatively, 
a count of the eigenvalues can be developed in a local neighborhood of the bifurcation curves, such as the ones shown 
on Fig. 01 These computations can be developed analytically, with the perturbation series expansions [islfl^ lr ^ . 



A. Scalar solitons 



For simplicity and without loss of generality, we set Ua = in our analytical computations. First, we study the 
stability of scalar solitons, when u = (f>{x, y) and t; = 0, where 4'{x, y) is a solution of the nonlinear eigenvalue problem: 

A0-/3(/. + y(a:,y)(,5 + 02)<^ = O, (6) 

We assume that there exists a ground state (positive definite) solution of the linear problem 

A</'o-/3o0o + '5V^(2;,2/)'/'o = O. 

with the propagation constant Pq. Applying the local bifurcation analysis for the nonlinear ground state, we look for 
the solutions in the asymptotic form, (f) = e[0o + e^'/'2 + 0(e^)] and /3 = /3o + e^/?2 + 0(e'*), and obtain the result 

"2- TT -7—\ > ^' ^'1 

(00, 00 ) 

where the inequality follows from the fact that V{x, y) is non-negative. Therefore, the soliton power (squared 
norm) p{(3) — 1|0||^2 = (0,0) is an increasing function of the propagation constant /3 near the bifurcation point 

/3 = /3o: 

f^-J^,Oi^,>. (S, 

Stability for the scalar solitons is determined by the linear eigenvalue problem, L^u = —Xw and L-w = Xu, where 
the linear operators L± are defined as, 

L+ = P-A-Vix,y){S + 3^''{x,y)) 
L_ = i3-A-V{x,y){8^^^(x,y)). 

If (j){x,y) > for all {x,y) € R^, then L- is non-negative with the zero eigenvalue L_0 = due to the gauge 
invariance. We shall consider the number of negative eigenvalues of and apply the earlier results for one-dimensional 
solitons 0. It is clear that L_|_ must have at least one negative eigenvalue, since 

(0,L+0)--2(02,y(a;,y)02)<o (9) 

When [3 = Pi) and = 0, the operator L-|_ has a simple zero eigenvalue and no negative eigenvalues. Therefore, 
according to the perturbation theory, the operator has exactly one negative eigenvalue for [3 > (3q near the local 
bifurcation threshold. The condition for applicability of the Vakhitov-Kolokolov criterion is satisfied and it suggests 
stability of scalar solitons at least near the bifurcation point. 

Numerically, we have checked that the number of negative eigenvalues of does not change and the slope of 
dP/d(3 is always positive, as shown in the example presented in Fig. ^c). Therefore, the scalar optical solitons in 
PCF is stable everywhere for (3 > (3q. 



B. Vector solitons 



Next, we study stability of vector solitons, when u = ^{x,y) and v ~ ^'(a;,?/), where ^{x,y) and ^'(a;,?/) are 
real-valued positive solutions of the coupled nonlinear eigenvalue problem: 



^$ = A$ + V{x, y) ((5 + $2 + 
7* = A* V{x, y) {6 + /i*^ + 



(10) 
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We consider a local bifurcation of the vector soliton from the scalar one, and look for solutions in the asymptotic 
form, $ = $0 + £^^"2 + 0(e'*), = e(5'o + £^^'2 + 0(6^^)), and also expand the eigenvalue, 7 = 7o + e^72 + 0(e''), where 
(3 is an arbitrary parameter, such that f3 > (3q. The function $0 = 'Pi^ill) satisfies the nonlinear eigenvalue problem 
® for a scalar soliton. Function 5*0 = ip{x, y) is a ground state solution of the linear eigenvalue problem: 

Lo^ = (70 - A - V{x, y) {5 + /i^')) = 0, (11) 

where 70 is a function of parameters ((S, /^) and the propagation constant [3. The problem for (^2{x,y), -L+$2 = 
^V{x,y)4>'ij}^ , is always solvable, in the assumption that the operator for the scalar soliton has one negative and 
no zero eigenvalues for any (3 > (3q. Finally, from the solvability condition of the linear inhomogeneous problem, 

Lq^2 - 2fiVix, y)(l}ip^2 + Vix, y)^^ - 72V', (12) 

we derive that 

2fi{^p^,V{x,y)cj,<^2) + {i^^V{x,y)ij^) 



72 



Numerical results show that 72 > for /i < 1 (i.e. the bifurcation occurs from the lower boundary of the existence 
domain on the plane (/3,7)), 72 < for fi > 1 (i.e. bifurcation occurs from the upper boundary of the existence 
domain), and 72 = for /i = 1 (i.e. the existence domain shrinks on the diagonal 7 = /3 > /3o and 70 = (3q.) 

We compute the Hessian matrix of derivatives of individual powers P = ($, $) and Q = (vp, ^) with respect to 
parameters P and 7. Let p = (0, </<) and assume that p'{(3) > for scalar soliton with (3 > Pq. Near the local 
bifurcation threshold at 7 = 70, we have: 

|e = ?MH),0(.') = §. |3 = <M)+0(.»,, (13) 

07 72 dp 72 

such that the determinant of the Hessian matrix is 



^(/3,7)=(^^^«^-l^+0(e^). (14) 



\2 

2 

72 I2 



When /i < 1, we have 72 > and the determinant may change the sign. Numerical results show that D > for 
Po < P < P* and D < for /? > /3, near the local bifurcation boundary 7 = 79. 

When /i > 1, we have 72 < and the determinant is always negative near the local bifurcation boundary 7 = 70- 
When yit = 1, we have P = j and P = Q, such that D = 0. 

With the standard linearization, the stability problem for vector solitons reduces to the matrix eigenvalue problem 
L+u = —Aw and L-w = Au, where u is a two-vector of real parts of the perturbation and w is a two-vector of 
imaginary parts of the perturbation, for a real eigenvalue A. The matrix Schrodinger operators are 

+ ~l -2^V"$* 7- A-\/((5-H3«'2-K^$2) 



f _ / /3 - A - y((5 $2 + ^^-2) \ 

\ j-A-ViS + ^^ + fi^^) J 

Since L_ is a diagonal composition of two scalar Schrodinger operators, each has a simple zero eigenvalue with the 
ground state $ and therefore, the operator i_ is non-negative. Therefore, stability of fundamental vector solitons 
is determined by the number of negative eigenvalues of the matrix operator Z+, similar to Ref. [l^ . 

We compute the number of negative eigenvalues of the operator L+ near the local bifurcation point. When 7 = 70 
and 5* = 0, we have 

such that the operator ij+ has exactly one negative eigenvalue (by the assumption that the scalar soliton is stable 
for P > Po) and the operator Lq has a simple zero eigenvalue with the eigenfunction ip. We study bifurcation of the 
simple zero eigenvalue of Lq for 7 7^ 70. Using the same small parameter e as in the local bifurcation analysis, we 
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are looking for solution of the eigenvalue problem Z+u = Au by the regular perturbation theory: ui = eUi + 0(e'^), 
U2=tp + + 0(e''), and A = 6^X2 + 0(e4). 

By algorithmic computations of the regular perturbation theory, we have the linear inhomogeneous problem for the 
first-order correction, L+Ui — 2iJ,V{x,y)(j)ip'^ , which is solvable with the solution Ui = 2^2- Furthermore, we have 
the linear inhomogeneous problem for U2, 

L0U2 = (A2 - 72)^ + 2^lV{x, y)H^2 + 3F(x, y)^^ + 2^1/(a;, y)HUi, 

with the solvability condition: 

6M(V'^ V{x, z/)0$2) + 3(V>^ V{x, y)i>^) 

^'^'^ W¥) ^ 

When /i < 1, we have 72 > 0, such that the zero eigenvalue of becomes a negative eigenvalue of L^. As a 
result, we have two negative eigenvalues of i+ near the local bifurcation boundary. Since p'{(3) > 0, we have two 
positive eigenvalues of the Hessian matrix when D > and one positive eigenvalue when D < Q. In the former 
case, the vector soliton is stable, while it is unstable in the latter case. Therefore, the boundary of _D = separate 
the domains of stability and instability of vector solitons on the plane (/3, 7) in the assumption that the number of 
negative eigenvalues of remains unchanged in the entire existence domain. 

When /I > 1, we have 72 < 0, such that the zero eigenvalue of Lq becomes a positive eigenvalue of L+. As a result, 
we have only one negative eigenvalue of L+ near the local bifurcation boundary. In the same region, we have exactly 
one positive eigenvalue of the Hessian matrix, since _D < 0. Therefore, the vector soliton is stable near the local 
bifurcation boundary. Numerics show that there exists a curve = in the existence domain (see Fig. ^Jb)), where 
the positive eigenvalue of the Hessian matrix crosses zero and becomes negative eigenvalue. These curve approaches 
the bifurcation curves asymptotically for large (/3,7), since Z? < on the bifurcation curves. In the assumption 
that the number of negative eigenvalues of remains unchanged in the entire existence domain, the curve _D = 
separates the stability and instability domains. 

When /i = 1, we have 7 = /3 and the zero eigenvalue of Lq is preserved as the zero eigenvalue of in the entire 
existence domain P — ^ > (3q. This additional eigenvalue is related to an arbitrary polarization of the vector soliton 
in the case ^ = 1: ^ — cos 6 (p and ^' = sin6' 0, where solves the scalar problem (0. The operator L+ always has 
a single negative eigenvalue (since we have verified numerically that L+ has a single negative eigenvalue for (3 > Pq). 
Therefore, the vector soliton must be linearly stable in the case /i = 1 for any f3 > (3o (excluding the limit (3 00). 



VI. CONCLUSIONS 



We have demonstrated that stable two-dimensional vector solitons can be supported by a nonlinear PCF structure 
with the Kerr nonlinearity. They constitute a class of two-component spatially localized modes that bifurcate from 
their one-component scalar counterparts and are described by two independent parameters. Both scalar and vector 
solitons provide a generalization of the guided mode trapped in the PCF core to the nonlinear case, being confined 
by both linear and self-induced nonlinear refractive indices. The periodic PCF environment provides also an effective 
stabilization mechanism for these localized modes, in a sharp contrast with an entirely homogeneous nonlinear Kerr 
medium where both scalar and vector spatial solitons are unstable and may undergo the collapse instability. We have 
applied the analytical matrix criterion for stability of these PCF vector solitons, and have verified that this criterion 
is confirmed by the direct simulations of the soliton dynamics. 



VII. ACKNOWLEDGMENTS 



This work has been supported in part by the Australian Research Council. JRS acknowledges a visiting fellowship 
granted by the Direccion Xeral de Investigacion e Desenvolvemento of Xunta de Galicia (Spain). Both JRS and DEP 
thank Nonlinear Physics Center at the Australian National University for a warm hospitality during their stay in 
Canberra. 



[1] P. Russell, Science 299, 358 (2003) 



9 



[2] J.C. Knight, Nature 424, 847 (2003). 

[3] Yu.S. Kivshar and G.P. Agrawal, Optical Sohtons: From Fibers to Photonic Crystals (Academic, San Diego, 2003). 
[4] Y. V. Kartashov, A. S. Zolcnina, V. A. Vysloukh, and L. Torncr, Phys. Rev. E 70, 066623 (2004). 
[5] Z. Chen, A. Bezryadina, I. Makasyuk, and J. Yang, Opt. Lett. 29, 1656 (2004). 
[6] S.F. Mingaleev, and Yu.S. Kivshar, Phys. Rev. Lett. 86, 5474 (2001). 

[7] A. Ferrando, M. Zacares, P. Fernandez de Cordoba, D. Binosi, and J. A. Monsoriu, Opt. Exp. 11, 452 (2003). 
[8] A. Ferrando, M. Zacares, P. Fernandez de Cordoba, D. Binosi, and J. A. Monsoriu, Opt. Exp. 12, 817 (2004). 
[9] J.N. Malmborg, A.H. Carlsson, D. Anderson, M. Lisak, E.A. Ostrovskaya, and Yu.S. Kivshar, Opt. Lett. 25, 643 (2000). 
[10] J. E. Dennis and R. B. Schnabcl, Numerical Methods for Unconstrained Optimization and Nonlinear Equations (Englcwood 
Chffs, NJ: Prcntice-Hall, 1983) 

[11] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical recipes in C (Cambridge University Press, 

Cambridge, 2nd ed. 1992). 
[12] C. C. Paige and M. A. Saunders, SIAM J. Numer. AnaL 12, 617 (1975). 
[13] D.E. Pelinovsky and Yu.S. Kivshar, Phys. Rev. E 62, 8668 (2000). 
[14] D.E. PcUnovsky, Proc. Roy. See. Lond. A 461, 783 (2005) 
[15] D.E. PeUnovsky and J. Yang, Stud. Appl. Math. 105, 245 (2000). 
[16] J. Yang and D.E. Pohnovsky, Phys. Rev. E 67, 016608 (2003). 
[17] D.E. Pehnovsky and J. Yang, Stud. Appl. Math., in print (2005). 



